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Abstract 

In this note we shall introduce a simple, effective numerical method for solving 
partial differential equations for scalar and vector-valued data defined on surfaces. 
Even though we shall follow the traditional way to approximate the regular surfaces 
under consideration by triangular meshes, the key idea of our algorithm is to develop 
an intrinsic and unified way to compute directly the partial derivatives of functions 
defined on triangular meshes. We shall present examples in computer graphics and 
. image processing applications. 
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1 Introduction 



Numerical approaches to solve partial differential equations (PDE's) on sur- 
faces have received growing interest over last decade. However, they are still 
not well-understood. Partial differential equations need to be solved intrinsi- 
cally and numerically for data defined on 3D surfaces in many applications. 
For instance, such examples exist in texture synthesis (Turk[12j, Witkin and 
Kass[13j), vector field visualization (Diewald, Preufer and Rumpf[?]), weather- 
ing (Dorsey and Hanrahan[6j) and cell-biology (Ayton, McWhirter, McMurty 
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and Voth 2005). Usually, surfaces are presented by triangular or polygonal 
forms. Partial differential equations are then solved on these triangular or 
polygonal meshes with data defined on them. The use of triangular or polygo- 
nal meshes is very popular in all areas dealing with 3D models. However, it has 
not yet been a widely accepted method to compute differential characteristics 
such as principal directions, curvatures and Laplacians (Chen and Wu[2j, Wu, 
Chen and Chi 03] ? Taubin[9,10]). This is because that there is no unified, 
simple and effective method to compute these first and second order differen- 
tial characteristics of the triangular or polygonal surface and to solve PDE's 
for data defined on triangular or polygonal meshes. In Chen, Chi and Wu[5], 
the authors proposed a new intrinsic simple algorithm to handle this difficulty. 
In this note, we shall use this new technique to solve PDE's on surfaces. 

In Bertalmio, Cheng, Osher and Sapiro[T] proposed a framework, the implicit 
surface algorithm, to solve variational problems and PDE's for scalar and 
vector-valued data defined on surf aces. Their key idea is to use, instead of a 
triangular or polygonal representation, an implicit representation. The surface 
under consideration is the zero-level set of a higher dimensional embedding 
function. Then they smoothly extend the original data on the surface to the 
3D domain, adapt the PDE's accordingly, and implement all the numerical 
computations on the fixed Cartesian grid corresponding to the embedding 
function. The advantage of their method is the use of the Cartesian grid instead 
of a triangular mesh for the numerical implementation. 



2 Our intrinsic algorithm for solving PDE's on surfaces 

In ths section we first propose our discrete intrinsic algorithm for solving 
PDE's on regular surfaces. We divide our algorithm into two main steps: First, 
we approximate the given surface by a suitable triangular mesh according to 
the accuracy demand. Second, we use our new intrinsic differential method 
developed in Chen, Chi and Wu[5j to compute the numerical PDE on the 
fixed triangular mesh. The first step is now easy to implement since one can 
find some good and efficient algorithms in the public domain. The difficult part 
lies in the second step. Namely, how can one effectively compute differential 
quantities on functions on a triangular mesh? 

Next, we shall compare our algorithm with the implicit algorithm proposed 
by Bertalmio, Cheng, Osher and Sapiro[I]. We list the key steps of these two 
algorithms about solving PDE's on surfaces as follows. 

One can tell from this comparison that our method is much simpler and more 
intrinsic. In many applications, one usually starts with triangular meshes in- 
stead of regular surfaces. In this case, we do not need Step 1 in our intrinsic 
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Our intrinsic algorithm 


The implicit surface algorithm 


1. Obtain a triangular mesh approx- 
imated to the given surface. 


1. Obtain an implicit representation 
of the given fixed surface. 




2. Extend smoothly the data on the 
surface to the 3D volume 




3. Adapt PDE's accordingly 


2. Use our new intrinsic differential 
method to compute the numerical 
PDE's on the fixed triangular mesh. 


4. Perform all the computations on 
the fixed Cartesian grid correspond- 
ing to the embedding function. 



algorithm. However, in the implicit surface algorithm, one will need one ex- 
tra processing step: Construct an accurate implicit surface from a triangular 
mesh. Note that the triangular mesh may compose of a lot of triangles. This 
will cost large computations to obtain the accurate implicit surface. 

Next, we shall describe a new, simple and effective method to define the dis- 
crete gradient and the discrete LB operator on functions on a triangular mesh. 
In order to do so, we first recall the gradient and the LB operator on functions 
in a regular surface E in the 3D Euclidean space 1R 3 . 



2. 1 Gradient and LB operators on regular surfaces 



We consider a parameterization x : U — > E at a point p, where U is an open 
subset of the 2D Euclidean space M 2 . We can choose, at each point q of x{U) ) 
a unit normal vector N(q). The map TV : x{U) — ► S 2 is the local Gauss map 
from an open subset of the regular surface E to the unit sphere in the 3D 
Euclidean space M 3 . The Gauss map TV is differentiate. Denote the tangent 
space of E at the point p by TE p = {v e K 3 : V-LN(p)}. The tangent space 
TYip is a linear space spanned by {x U) x v } where v are coordinates for U. 

The gradient Vg of a smooth function g on E can be computed from 



^ 9uG - g v F g v E-g u F 

EG-F* Xu ^ EG-F* Xv [ } 

where E ) F, and G are the coefficients of the first fundamental form and 
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dg(x(u,v)) 
du 

dg(x(u,v)) 
dv 



(2) 



3 



See do Carmo[7j for the details. Note that the gradient Vg assigns to each 
point q in E a tangent vector Vg{q) such that we have for all v G TE g , 

W9W>v) q = — - — |U=o (3) 
where the smooth curve 7(f) is in E with 7(0) and 7 / (0) = v. 
The LB operator A acting on the function g is defined by the integral duality 



(Ag,<l>) = -(Vg,V<f>) (4) 

for all smooth function (f) on E. A direct computation yields the following local 
representation for the LB operator on a smooth function g: 



Vg 



Veg—f' 2 



G 



d 



du V ^EG-F 2 du> du\ ^EG-F 2 dv ) 



+ I \A( E_ 



F 2 dv ) dv^VEG-F 2 du\ 



I) _ d_( 



.5a 



(5) 



To move from regular surfaces to triangular meshes, one need to avoid the 
problem of local parametrization x around a vertex p. In other word, one 
does not have the fist fundamental form F, G and their derivatives for the 
computation of the gradient and the Laplacian operator of a function on a 
triangular mesh. To handle this problem, we give a novel method in Chen, 
Chi and Wu[5j to compute these differential quantities. The primary ideas 
were developed in Chen, Chi and Wu[4j where we try to estimate the discrete 
partial derivatives for 2D scattered data points. 



2.2 A new discrete algorithm: local tangential lifting (LTL) method 



In this section we shall describe a unified, simple and effective method to 
define the discrete gradient and the discrete Laplacian operator on functions 
on a triangular mesh. The primary ideas were developed in Chen, Chi and 
Wu[3||l] where we try to estimate the discrete partial derivatives of functions 
on 2D scattered data points. Indeed, the method that we shall use to develop 
our algorithm is divided into two main steps: first we lift the 1-neighborhood 
points to the tangent space and obtain a local tangential polygon. Second, we 
use some geometric idea to lift functions and vectors to the tangent space and 
then we can compute their derivatives in the 2D tangent space. This means 
that the lifting process allows us to reduce the 2D curved surface problem to 
the 2D Euclidean problem and hence the methods in [4J and[8j can be applied. 

Consider a triangular mesh S = (V, F), where V — {vi\l < i < ny} is the list 
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Fig. 1. The tangential polygon P{v). 

vertices and F = {fk\l < k < rip} is the list of triangles. Next, we introduce 
the notion of the local tangential polygon P(v) at the vertex v of S as follows: 

(1) The normal vector N(v) at the vertex v in S is given by 

^2feT(v) VfNf 



N(v) = 



(6) 



where Nf is the unit normal to a triangle face / and the centroid weight 
is given in [2] by 



UJf 



\\G 



f- 



(7) 



f£T(v) \\G 



f 



Here, Gf is the centroid of the triangle face / determined by 



(2) The tangent plane TS(v) of 5 1 at i> is now determined by 

T5(v) = {w£ R 3 \w±N(v)}. 



(9) 



(3) The local tangential polygon P{v) of v in TS(v) is formed by the vertices 
5$ which is the lifting vertex of Vi adjacent to v in S. 



Vi = (vi — v)— < Vi — N(v) > N(v). 



(10) 



as in figure [TJ 



Let ft be a function on V. We will lift locally the function h to a function, 
denoted by h V) on the vertices in P{v) by simply setting 



h v (vi) = h(vi). 



fill 
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And h v (0) = h(v) where is the origin of TS(v). One can then extend the 
function h v to a piecewise linear function, still denoted by h V) on P(v) as 
follows. 

Consider a face / with vertices v, Vi and vj in F. We obtain a lifting face / 
with vertices 0, vi and Vj in P(v). Every point p in / can be written as a linear 
combination of Vi and Vj . That is, p = avi + bvj where a, b > and a + b < 1. 
Then we define 



h v (p) = ah v (vi) + bh v (vj) + (1 - a - b)h v (6). 



(12) 



Hence, the extended function h v is affine on each triangle / of P{v) and is 
differentiate on /. The gradient V(h v )j of h v at the origin can be obtained 
by 

V(h v )j(6) = avi + (3vj (13) 
where the coefficients a and f3 satisfy the relations: 



h(vi) - h v (6) =< (Vh v ) f (6),Vi > 
^ h(vj) - h v (6) =< (V^) 7 (0),^ > 



(14) 



As easy computation gives 



a 



\0, 



< Vi,Vi > < Vi,Vj > 
^< Vi,Vj > < Vj,Vj > i 



h v (vi) - h v (6) 



(15) 



To obtain the gradient Vh(v) of h on S at the vertex v, we use again the 
weighted combination method. Namely, we set 

Vh(v) = (Vh v )(6)= J2 ^f(yh v ) f i0) (16) 

feP(v) 

with 



\\G~W 2 

"f = v f i (17) 
where Gj is the centroid of the lifting triangle face / and is determined by 

G,= H^. (18) 



Next we explain how to obtain a good discrete Laplacian Ah(y) of a function h 
on the triangular mesh S. From the discussions above, we obtain the gradient 
Vh(v) of h on S at each vertex v. We can use the method of parallel transport 
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Fig. 2. parallel transport 

to lift the vector Vh{vi) at V{ to a vector Vh(vi) in the tangential space TS(v). 
The idea is to define a orthonormal linear map from TSivi) to TS(v). To do 
so, we choose an orthonormal basis for TS(v) by 



N(v) 



_ (vi—v) — <Vi—v,N(v)>N(v) 
~~ \\(vi-v)-<Vi-v,N(v)>N(v)\\ 



(19) 



e 2 = N(v) x e 1 



The corresponding orthonormal basis for TS(vi) is then given by 



N(vi 



ei 



__ (t>»-t))-<^-t>,Affa)>Af(t>i) 
\\(vi-v)-<Vi-v,N(vi)>N(vi)\\ 



(20) 



e 2 = iV(^) x ei 
Then, the linear map L of the parallel transport is given by 

L(w) = ae 1 + be 2 G TS(v) 
for u> = aei + be 2 in TS(vi). See figure El 
In this way, we can set the tangential gradient Vh at Vi by 

V^) = L(Vh( Vi )). 



(21) 



(22) 



Hence we obtain a tangential gradient Vh of /i at each vertex Vi in the tan- 
gential polygon P(v) and we also set 



V/i(0) = V/i(u). 



(23) 



See figure [31 



7 



Fig. 3. Vh 

Fix an orthonormal basis {ei, e 2 } for TS(v). The tangential gradient Vh can 
be written as 

Vh(vi) = a(vi)ei + b(v z )e 2 . (24) 

The coefficients a(^) and b(vi) can now be viewed as functions on the vertices 
Vi of the tangential polygon P{v). As before, we can then obtain their gradients 
(Va)(0) and (V6)(0) at origin. Namely, 




(Va)(0) = E/ e p ( ,) ^/(Va^O) 



(V6)(0) = E/ 6 p W w/(V&)/(0) 
with the centroid weights a;^ as in (IT7I) . 

Put them in the matrix form to give 



(Va)(0) = a n ei + a 2 ie 2 

( 26 ) 

(V6)(0) = a 12 e 1 + a 22 e 2 . 



Therefore we have the Laplacian Ah(v) of h at the vertex v of S 1 : 

A/i(u) = an + a 22 . (27) 



Theoretically the definition of the Laplacian Ah(v) is independent of the 
choice of the orthonormal basis {ei, e 2 }. 
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3 Linear diffusion 



As a simple example to illustrate our new algorithm, let us consider a linear 
diffusion equation on a regular surface E: 

u t - Au = g on Ex/ (28) 
for u : S x / ^ R, / C K, where A is the surface Laplacian on E and 
g:Ex/^Misa smooth function on E. For the numerical implementation 
of our intrinsic algorithm, we take the regular surface E to be (1) the unit 
sphere S 2 or (2) a torus T 2 . In the case of the sphere 5 2 , we consider the 
function 

g(x) = X\ for x = (xi, x 2 , x 3 ) G S 2 . (29) 

Figure @] and figure [5] give the solution of (1251) and (|29l) with initial functions 
?x(x, 0) = 0. Different time steps are shown until the stationary solution is 
reached. 

Consider the torus T 2 a r ^ = ((a+r cosx) cos?/, (a+r cosx) sin?/, sinx) for x, ?/ G 
[0, 27r] with a > r > 0. we take a = 2, r = 1 and choose the function 

g(x,y)=x for x, ?/ G [0, 27r] (30) 

Figure [6] and figure [7| give the solution of (1251) and (|30l) with initial functions 
?z(x, 0) = 0. As above, different timesteps are depicted until the stationary 
solution is reached. 



4 Reaction-diffusion textures 

The original idea about how reaction-diffusion equations can be used to cre- 
ate patterns was first introduced in (TuringfTTj). The basic idea is to have a 
number of chemicals that diffuse at different rates and that react with each 
others. After the works of (Turk[T2]. Witkin and Kass[13]), the use of reaction- 
diffusion equations for texture synthesis attracted a lot of attentions in com- 
puter graphics. Turk , Witkin and Kass used these equations for planar tex- 
tures and textures on surfaces. Then the patterns are analyzed by assigning a 
brightness value to the concentration of one of the " chemicals" . 

Consider two chemicals U\ and u 2 on a surface E. In a simple isotropic model, 
we have 
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< 



dui 
dt 



f(u u u 2 ) + ttAli! 
g{u 1 ,u 2 ) + (3Au 2 



(31) 



du2 - 
dt ~ 



where a and (3 are two constants representing the diffusion rates and / and g 
are functions that describe the reaction. For simple isotropic patterns, Turing 
chose the functions / and g to be 

f(u u u 2 ) = s(16 - u ± u 2 ) 
< (62) 

< g{u X) u 2 ) = s{uiy 2 - y 2 - 7) 

where s is a constant and 7 is a random function giving the irregularities in 
the chemical concentration. 

By using our intrinsic method described in the previous section, we can easily 
generate textures on surfaces without the elaborated schemes employed in 
(Turk 1991, Witkin and Kass 1991). 

In the case of the sphere 5 2 , figure [8] and figure [9] give the solution of (I3TI) 
and (|32l) with initial functions and Different timesteps are shown until the 
stationary solution is reached. 

On the torus Tj? r \ = ((a+r cos x) cos y, (a+r cos x) sin y, sin x) for x, y e [0, 2tt] 
figure [10] and figure [TT] give the solution of (I3TI) and (1321) with initial functions 
u\{x, 0) ^2(^5 0) = 1 and a = l,/? = s = 2, 7 = 0. Different timesteps are 
shown until the stationary solution is reached. 
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Fig. 10. torus 




